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This paper considers the stability of liquid metal drops subject to a high-frequency AC 
magnetic field. An energy variation principle is derived in terms of the surface integral 
of the scalar magnetic potential. This principle is applied to a thin perfectly conducting 
liquid disk, which is used to model the drops constrained in a horizontal gap between 
two parallel insulating plates. Firstly, the stability of a circular disk is analysed with 
respect to small-amplitude harmonic edge perturbations. Analytical solution shows that 
the edge deformations with the azimuthal wavenumbers m = 2, 3, 4, . . . start to develop 
as the magnetic Bond number exceeds the critical threshold Bm c = 3ir(m + l)/2. The 
most unstable is m = 2 mode, which corresponds to an elliptical deformation. Secondly, 
strongly deformed equilibrium shapes are modelled numerically by minimising the as- 
sociated energy in combination with the solution of a surface integral equation for the 
scalar magnetic potential on an unstructured triangular mesh. The edge instability is 
found to result in the equilibrium shapes of either two- or three-fold rotational symmetry 
depending on the magnetic field strength and the initial perturbation. The shapes of 
higher rotational symmetries are unstable and fall back to one of these two basic states. 
The developed method is both efficient and accurate enough for modelling of strongly 
deformed drop shapes. 



1. Introduction 

In several metallurgical processes such as, for example, the levitation melting and cold- 
crucible, where the induction heating is used, the surface of liquid metal is subject to 
AC magnetic field. In such a way, the metal can be not only molten but also evaporated 
provided that the heating power is high enough (Baptiste et al. 2007). Induction heating 
is accompanied with a pinch effect, which can significantly deform the surface of liquid 
metal. When a sufficiently strong magnetic field is applied, surface sometimes becomes 
asymmetric and even strongly irregular (Fautrelle et al. 2007). This phenomenon is of 
primary importance for the induction heating of liquid metals because it may have an 
adverse effect on the heating efficiency and eventually limit the power density the liquid 
metal can dissipate. Such a surface instability has been observed first by Perrier et al. 
(2003) on a circular layer of Gallium in a mid-frequency AC magnetic field. Analogous 
instability was studied also by Mohring et al. (2005) on the free surface of InGaSn 
melt in the annulus placed under a ring-like circular coil and fed by an alternating 
current with the frequency in the range of 20 — 50 kHz. As the current amplitude exceeds 
a certain critical value, which depends mainly on the annulus width, an initially flat 
surface acquires a static wavy deformation. At a higher critical current, the deformation 
rapidly increases and becomes unsteady. In contrast to Perrier et al. (2003), who observe 
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only static deformations, Kocourek et al. (2006) find a circular sessile drop of InGaSn 
melt first to squeeze radially with various shape oscillations to set in as the strength of 
a 20 kHz AC magnetic field is gradually increased. In the former experiment, the surface 
of liquid metal was exposed to the air and, thus, heavily oxidised that constrained its 
motion. In the latter experiment, oxidation was prevented by covering the drop by a 
diluted HC1 solution. Later on Conrath et al. (2006); Conrath (2007) found static 
shape deformations when the drop was constrained in a horizontal gap between two 
parallel plates. Irregular static surface shapes have been observed also by Hinaje et al. 
(2006a) on a layer of PbSn alloy covering the bottom of a cylindrical container. The metal 
layer, which was constrained by the lateral walls of the container and heavily oxidised at 
the top, broke up revealing the bottom of the container as the strength of a 4 Khz AC 
magnetic field exceeded a certain critical value. The authors also attempted to model this 
process numerically using a surface integral equation derived from Green's third identity. 
This approach, however, is not applicable to thin sheets, for which the double layer 
contribution vanishes. In a subsequent paper, Hinaje et al. (2006b) devised a simplified 
electrotechnical model, which provided a rough estimate of equilibrium shapes. A simple 
theoretical model for this type of instability was introduced by Priede et al. (2006) , who 
analysed the linear stability of the edge of liquid metal layer, which was treated as a 
perfectly conducting thin liquid sheet in a transverse AC magnetic field. This allowed 
the authors to determine the wavenumber of the most dangerous perturbation and the 
critical field strength at which the instability develops in a reasonable agreement with 
the observations of Mohring et al. (2005). 

In this paper, an energy variation principle is derived for the equilibrium shapes that 
develop from the edge pinch instability of flat liquid metal drops, which are modelled as 
thin perfectly conducting liquid sheets. Firstly, the stability of a circular disk is analysed 
with respect to small-amplitude harmonic edge perturbations. Analytical solution shows 
that the edge deformations with the azimuthal wavenumbers m = 2,3,4,... start to 
grow as the magnetic Bond number exceeds the critical threshold Bm c = 37r(m + I)/2. 
The most unstable is m = 2 mode, which corresponds to an elliptical deformation. 
Secondly, strongly deformed equilibrium shapes are modelled numerically by minimising 
the associated energy. The electromagnetic problem is formulated in terms of the surface 
integral equation for the scalar magnetic potential, which is solved numerically on an 
unstructured triangular mesh covering the surface of the drop. The edge instability is 
found to result in the equilibrium shapes of either two- (m = 2) or three-fold (m = 
3) rotational symmetry depending on the initial perturbation and the magnetic field 
strength. Although the associated energy of m = 3 shapes is higher than that of m = 2 
ones at the same magnetic field strength, both shapes are separated by a positive energy 
barrier. This, however, is not the case for equilibrium shapes of higher order symmetries. 
Although these shapes can be obtained numerically, they turn out to be unstable with 
respect to small amplitude perturbations of two- or threefold rotational symmetries, 
which make them fall back to one of the two basic states. 

This paper is organised as follows. In §2, the problem is formulated and the energy 
variation principle derived in terms of the integral of the scalar magnetic potential over 
the drop surface. This principle is applied in §3 to obtain an analytical solution for the 
stability of a circular disk with respect to small-amplitude harmonic edge perturbations. 
Specific mathematical details of the solution are given in Appendix A. In §4, numerical 
method is described and validated against the previous analytical solution. Numerical 
results are presented in §4. The paper is concluded with a summary and discussion in §5. 




2. Formulation of problem 

Consider a drop of liquid metal with the characteristic size Rq, electrical conductivity 
a, surface tension 7 and density p submitted to an AC magnetic field with the spatial 
amplitude distribution B(r), , as shown in figure f. The AC frequency u> is assumed 
so high that the penetration depth of the magnetic field into the drop 5 <~ (fiQauj)^ 1 ^ 2 , 
where p is the vacuum permeability, is negligible with respect to Rq- In this perfect 
conductor approximation, the magnetic field is tangential to the drop surface S, 

B n \ s = 0, (2.1) 

and the electromagnetic force effectively acts on the surface as the time-averaged mag- 
netic pressure, 

B 2 

Pm = ~. ■ 

Equilibrium shape of the drop is determined by the normal stress balance 

Ph -p c -Pm\ s = 0, (2.2) 



where — pgr and p c = 7 V n are the hydrostatic and capillary pressures, respect 
n is the outward surface normal and g is the gravitational acceleration. Multiplying (2.2) 
by n£, where £(r) is a virtual displacement field conserving the volume, and integrating 
over S, we obtain 

W g + W s + W m = 0, (2.3) 

where W g = p J s (g ■ r)£ • ds, W s = -7 J s {V ■ n)£ • ds and W m = J s B 2 £ ■ ds are the 
virtual works done by the gravitational, surface tension and magnetic forces, respectively. 
Since all of these forces, including the magnetic one in the perfect conductor approxima- 
tion, are conservative, the corresponding works can be expressed as the variations of the 
associated potential energies. Using the divergence theorem to change from the surface 
to volume integrals, and taking into account the incompressibility constraint V • £ = 
as well as the Lagrangian variation £ • V/ = Sf, we obtain W g = —SE g , W s = —SE S , 
and W m = —5E m , where 



E g = - f pg-rdV, (2.4) 
Jv 

E s = 7 5, (2.5) 

1 r -,2, 



E m - -— / B A dV (2.6) 

4^0 
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are the associated potential energies. The minus sign at the last integral is due to the 
integration over the outer volume V. For an equilibrium shape, (2.3) implies SE = 0, 
where 

E = E g + E s + E m (2.7) 

is the total associated energy. This derivation of the energy variation principle appears 
more straightforward than the original one by Sneyd & Moffatt (1982). Equilibrium 
shape of the drop corresponds to a stationary point of E, which, as usual, has to be a 
minimum for the equilibrium to be stable (Chandrasekhar 1961). 

2.1. Magnetic energy in a homogeneous external field 

Further, the external magnetic field B e is assumed homogeneous, which supposes the 
drop to be small compared to the inductor generating the field. This allows us to express 
the magnetic energy (2.6) by an integral over the drop surface as follows. The total 
magnetic field is a superposition of the external and induced fields B = B e + Bi. Outside 
the drop, we have Bi = -^V'l'i, where is the scalar potential of the induced magnetic 
field. Then (2.6) can be represented as 

E m = Eq + E\, (2.8) 

where 



E 1 = --^- [ B BidV = -( ViBe-ds, 
4 Mo Jv 4 Js^ 



(2.9) 



and Soo is a remote surface enclosing the drop at r — > oo. The part of the integral over 
the drop surface S vanishes because of the boundary condition (2.1). Since the induced 
magnetic field is supposed to fall off at large distances r — > oo as the dipole field with 
<~ 1/r 2 , the last integral converges to a non-zero value. The other contribution to the 
magnetic energy is 

E = --^- [ B-B e dV=--±- [ B\dV-±- [ Bi ■ B e dV, (2.10) 
4^o Jv 4/xo J v 4fi J v 

where the first integral represents the energy of the external magnetic field, which is 
constant, and thus negligible, however, formally it is infinite. Therefore, retaining only 
the second term in (2.10), we obtain 

E ° = \j~ v ' ( - Be * j) dv = ~\ S ^ iBe ■ ds + Eu (2 - n) 

where the integral is taken over the drop surface S with the outward normal direction. 
Now it remains to evaluate the integral in (2.9), which is determined by the dipole 
component of the induced field 

1 m ■ r 

**< P > = ^ — ' 

where m = i J s r x J ds is the dipole moment of the drop and J is the surface current 
density. The latter is related to the magnetic field by Ampere's integral current law, 
which applied to a small surface element results in 

J = — n x B\ s = V* x n\ s , (2.12) 
Mo 

where >I> = vl'e+^'i is the full scalar magnetic potential including also that of homogeneous 
external field 

* e = -— r B e . (2.13) 
Mo 
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Substituting these expressions into (2.9), after some algebra we obtain 

E 1 = ±B e -m=-^J *B e -da. (2.14) 
The integral above can be represented as 

£l = l(/ /(B «. d ,-M), 

where the second term, which is related to the energy of homogeneous external magnetic 
inside the drop of fixed volume, is constant and, thus, negligible again. By the same 
argument, "J/j in (2.11) can be substituted by Then the magnetic energy (2.8), apart 
from a constant contribution of the external field, can be written in terms of as 

E m = -\ [ VB e ■ ds + 2E l = -E u (2.15) 
4 Js 

where Ex is given by (2.14). 

2.2. Scalar magnetic potential 

There are two alternatives how to find the magnetic potential. First, the solenoidality 
constraint V • B = for a free-space magnetic field B = —^V^ results in 

V 2 * = 0, (2.16) 

which together with the boundary conditions, 

d„*| s = 0, and *U oc ^* e = -Mo 1 '"S e (2-17) 

governs ^ outside the drop. This formulation is used in § 3 for analytical treatment 
of small amplitude deformations of a circular disk by using a singular Taylor-series- 
type expansion around the basic state. Second, for efficient numerical solution, instead of 
(2.16), which has to be solved in the whole space outside the drop, it is more advantageous 
to use Biot-Savart law 

B(r) = B »ilRF xJ(r ' )dV ' (2 ' 18) 

where the prime denotes the integration point. Then the boundary condition (2.1) applied 
to (2.18) results in the surface integral equation defining $onS 

r' 

— • V*(r')dV = -n-B e , (2.19) 

which has to be solved for a given shape of the drop to obtain the surface distribution of 
<]/, which, in turn, defines the magnetic energy (2.14). Then equilibrium shape is found 
by minimising the total associated energy (2.7). 

2.3. Thin-drop model 

In the following, we focus on the case of a thin drop confined in a horizontal gap between 
to parallel insulating plates, as shown in figure 2. The drop is modelled by a perfectly 
conducting liquid sheet with the virtual displacements constrained to the plane of the 
sheet. The external magnetic field B e is perpendicular to the sheet. The surface enclosing 
the sheet consists of the top and bottom parts S+ and S- , both with the same area So but 
opposite normals. Taking into account also that the potential of the induced field changes 
the sign discontinuously by crossing the sheet, contributions from both surfaces in (2.14) 
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Figure 2. A model of liquid layer confined in the gap between two horizontal plates in a 

transverse AC magnetic field. 



are the same. This results in a factor of two at the front of the integrals in (2.14) and (2.19) 
when the integration is carried out only over the upper part of the sheet. Subsequently, 
we identify So and \& with the top surface and the potential at that surface, respectively. 
Note that the contribution of the transverse homogeneous field, whose potential (2.13) is 
constant along the sheet, vanishes in (2.14). For the layer of fixed thickness, gravitational 
energy is constant and, consequently, irrelevant in the variation of the total energy. Due 
to the volume conservation and fixed thickness, the horizontal area So is fixed, too. Then 
the variation of surface is caused only by the stretching of the perimeter P = § L dl, 
which determines the effective edge area S e — Plo and the corresponding surface energy 
E s — -)S e , where the arclength l over the edge is assumed to be fixed similar to the layer 
thickness itself (see figure 2). 

Subsequently, all variables are non-dimensionalised by choosing Ro, Bo, RoBo and 
7/0-R0 as the length, magnetic field, potential and energy scales, respectively. Then the 
dimensionless associated energy, which comprises a capillary contribution of the edge and 
the magnetic energy, can be written as 

E = I Al - -Bm / *ds, (2.20) 
Jl 3 J s 

where Bm = B^Rq/ (2/io7^o) is the magnetic Bond number based on the amplitude of AC 
magnetic field. Note that there is no difference between the induced and full magnetic 
field potentials in (2.20) when that of homogeneous field (2.13) is set to be zero along 
the sheet by a proper choice of additive constant. Actually, this difference is irrelevant 
because the contribution of homogeneous field in (2.13), as discussed in the previous 
section, is constant for incompressible liquid. For a flat sheet, (2.19) takes the following 
dimensionless form: 

^•V*(r')dV = -l. (2.21) 

For no electric current (2.12) to cross the edge L, r x n ■ J\ L = d T ^\ L = is required, 
which implies \I>| L =const, where r and r x n are the tangent and normal vectors to the 
edge, respectively, and n is the normal vector to the sheet. As discussed above, we can 
set 

*Il = «, (2.22) 
which ensures a zero potential for the homogeneous external field in the plane of the 
sheet. 
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3. Analytical solution for the stability of circular disk 

Here the approach developed above will be applied to analyse the stability of a circular 
liquid disk with the radius R = 1 + Ri + R 2 + . . . , where R x — Ri cos(m^) is a small 
perturbation with the amplitude /i^and the azimuthal wavenumber to, and i? 2 is a higher- 
order small correction to be determined later on. The potential is sought as = + 
*i + *2 + • • • , where 

*o(r?,0 = --7?[l + £arctan(0], (3.1) 

is the potential of circular disk presented in the angular and radial oblate spheroidal 
coordinates, < n < 1 and < £ < oo (Li et al. 2002), which are related with the 
cylindrical coordinates by 

z = r]£. 



Note that £ = corresponds to the plane of the disk z = 0, where r = y/l — vj 2 with 
n = corresponding to the edge of a circular disk at r = 1. The first-order perturbation 
of the potential vanishing away from the disk and satisfying the edge condition (2.22), 
which takes the form 



*i r _vi - -Ri 



= 4^ _1 U. ( 3 - 2 ) 



7T 



^(v,0 = -z(tT^ (3-4) 



dr 

can be written as 

* 1 (r) = R 1 i>?(r ] ,0, (3.3) 

where 

l-yf 

, i + e j v 2 +e 

The details of the solution above, which apart from slightly different notations are similar 
to those in Priede et al. 2006, can be found in Appendix A. Since the energy variation 
about the equilibrium state is expected to be quadratic in Ri, we need to consider 
also the next-order radius perturbation R 2 , which results from the area conservation 
S = J 2v /* r drd</> = tt(1 + Rl/2 + 2R 2 + . . .) as 

R 2 = -i??/4. (3.5) 

The second-order potential perturbation, for which the edge condition (2.22) takes the 
form 



*2 L-n = -ti\ — R2 



dr dr 2 dr 2 



2 

-_R 2 (1+cos(2to</>)) (mi]^ 1 + n~ 3 ) 

7T 



?7— >0 ' 



can be written as ^ 2 (r], £) = i? 2 (^(^O + ^"(^O cos(2to<^ . Subsequently, we will 
need only the first term of this expression 

which satisfies (Al) with to = 0, where *i(to£) is defined by (3.4). The second term 
above is obtained similarly to the first one by applying <9 2 to (3.1), as described in the 
last paragraph of Appendix A, which yields 

fro,,,)- 2 ^-^ + 3(1-^))] 
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At the disk surface, we have 



*o(r) = V 7 ! -r 2 

2 f m 

*y(r) = — 



tt°(r) 



7T Vl - r 2 ' 
2m+(l-r 2 )- 1 
7T Vl - r 2 



(3.7) 
(3.8) 

(3.9) 



It is important to note that (3.8) and (3.9) are singular at r = 1, which is the edge of the 
unperturbed disk. At the same time, the edge condition (2.22) implies the potential to 
be regular (zero) at the actual edge of the deformed disk. This implies that the solution 
above can be regularised by representing it in the radial coordinate f stretched with the 
radius of the deformed disk. Using the substitution 



r = Rf=(l + R)r, 



(3.10) 



where R = i?i + R 2 + . . . is the radius perturbation, and expanding the solution in power 
series of R up to the second order in R\, we obtain a solution of the same asymptotic 
accuracy, which is free of edge singularities 

- (Rf) 2 d 2 * 

*(r, 4>) = *(f (1 + R), 4>) « *(f , </>) + Rf— + ^-^—q^ + ■■■ = *(f , 0), 

where *(r, 0) = * M + (r) + R 2 (j> 2 {r) + *| m (r) cos(2m0)) + ... and 



*5"(r) = - 
tt°(r) = - 



2 r r ' 



3' 



7r \/T~— 

(m-l)(l-2r m )-r 2 2(r m - 1) 



(l_ r 2)3/2_ 

Then the magnetic energy term in (2.20) can be evaluated up the first order in R 2 as 

/ / *(r,^)rdrd</>- / R 2 *(f, 0)f df d<£ « -E m - R 2 E m . 
Jo Jo Jo Jo 

where E m and R 2 E m are the magnetic energies of circular disk and its leading-order 
perturbation defined by 



E„ 



-2tt / * (r)rdr = -, 



1: 

Jo 



*§(r) - 4*f(r) rdr = 4(m-l). 



£? TO = — 2n 

The surface energy term in (2.20) is evaluated as 

dZ w 2vr(l - i? 2 (m 2 - 1)) = E s + R 2 E S . 



(3.11) 
(3.12) 



Then the total energy variation is 
1 



SE 



-R 2 (E S + E m ) = -R 2 (m - 1) 7r(m + 1) 



-Bm 



Note that there is no energy variation for m = 1 mode, which corresponds to the shift of 
the disk as whole. Circular disk is stable with respect to small perturbation with m > 1 
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as long as its energy is at minimum, i.e. SE > 0. Since according to (3.5) R2 < 0, the 
stability condition for m = 2, 3, is satisfied as long as 

Bm< 3(m + l)7r/2. (3.13) 

The first unstable mode with m = 2, which corresponds to an elliptical deformation, 
appears as Bm exceeds the critical value 

Bm c - (3.14) 

This critical value is by a factor of 3 greater than the one found by our previous linear 
stability analysis (Priede et al. 2006). The cause of this discrepancy is discussed in the 
conclusion of the paper. 



4. Numerical solution 

This section introduces the numerical method which will be used subsequently to find 
equilibrium shapes of thin drops by the approach described in §2. Numerical solution 
will also be verified against the analytical results obtained in the previous section. To 
solve (2.21) with the edge condition (2.22), which define \I/ over the drop surface S, the 
latter is tiled into triangular elements as shown in figure 3(a). Triangulation is carried 
out as follows. Firstly we take a regular hexagon inscribed in the unit circle and tile it 
using equilateral triangles with the side length 1/N. Secondly, the hexagon is stretched 
radially to fit the unit circle. Then six points are discarded from the perimeter and 
the remaining 6(N — 1) points are redistributed uniformly against the midpoints of the 
previous radial level. This produces a more regular triangulation at the edge, which 
yields a slightly higher numerical accuracy. As a result, we obtain a triangular mesh 
with 6N 2 — 6 elements and 37V x (N + 1) — 5 vertices. Following the finite element 
approach, * is sought at the vertices and interpolated linearly within the elements. To 
determine '5 at the vertices, we need a corresponding number of equations, which are 
obtained by numerically approximating (2.21) at the inner points and applying the edge 
condition (2.22) at the peripheral points. The integral in (2.21) is represented as a sum of 
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10 15 20 25 30 35 40 

Radial number of elements, N 

Figure 4. The relative error in the magnetic energy of circular disk (3.11) against the radial 
number of elements N for four different Gaussian quadratures: (1) linear quadrature using only 
the centre point with the weight factor 1; (2a) and (2b) quadratic quadratures using three 
symmetric points with the barycentric coordinates (2/3, 1/6, 1/6) and (1/2, 1/2, 0), respectively, 
and the weight factors 1/3; (3) four-point cubic quadrature using the centre point with the weight 
factor —27/48 and three symmetric points with the barycentric coordinates (3/5, 1/5, 1/5) and 
weight factors 25/48 (Cowper 1973). 



integrals over separate elements, which are approximated by the Gaussian quadratures 
for triangles. Thus, for a given mesh n = (x i7 yi) 7 we obtain a system of linear equations 
with a dense matrix for unknown \I>j = \&(rj), which are found by the LU decomposition 
method. 

The convergence of the magnetic energy for circular disk, E m defined by (3.11), is 
shown in figure 4 against the radial number of elements N for four different quadratures. 
Accuracy is lower for the Gaussian quadratures with the evaluation points located closer 
to the mesh points. This is because of the integrand singularities encountered when the 
observation point belongs to the element over which the integral is evaluated. Subse- 
quently, we use a quadratic quadrature with three evaluation points located at the side 
midpoints of the element (curve 2b in figure 4), which provides the highest accuracy. For 
the linear elements used here, the integral in (2.21) can, in principle, be evaluated exactly. 
However, such an approach is not applicable because the singularities between adjacent 
elements do not cancel out when the current distribution is a piece- wise constant, as 
in this case, rather than continuous. Nevertheless, Gaussian quadratures still provide a 
reasonably accurate result also in this case. 

To verify the analytical solution obtained in the previous section, we first restrict the 
drop shape to an ellipse denned parametrically by the mesh point coordinates 

r i = {x° i R x ,y° i /R x ), (4.1) 

where (x°,y°) = r" are the mesh points for a circular disk and R x is a parameter 
defining the x-radius of ellipse. For a given Bm, equilibrium shape is found by using a 
Powell-type algorithm (Press et al. 1996) to minimise the associated energy (2.20) with 
respect to R x . For each R x , firstly, Vl/j is found by solving the system of linear equations 
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Magnetic Bond number, Bm 

Figure 5. Major radius 7? max = ma,x(R x , 1/R X ) versus the magnetic Bond number Bm for both 
stretched (R x > 1) and squeezed (7?^ < 1) along the x-axis ellipses with N = 16, 24 elements in 
the radial direction. The vertical dashed line shows the theoretical threshold value (3.14). 




for the corresponding distribution of mesh points (4.1). Secondly, integrals in (2.20) are 
evaluated numerically for the given distributions of Ti and ^i. 

As seen in figure 5, which shows the major radius of ellipse versus Bm, the critical value 
of Bm, by exceeding which the drop starts to deform, is slightly above its theoretical 
value (3.14). For N = 16 elements in the radial direction, the major radius slightly varies 
depending on whether the ellipse is stretched (R x > 1) or squeezed (R x < 1) along the 
x-axis. Although these two cases differ only by the orientation of the major axis of ellipse 
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Figure 7. Equilibrium shapes at Bm = 15 found with one azimuthal mode (M = 1) and N = 16 
elements in radial direction for the disks squeezed along the y- and :r-axes, which correspond to 
Rl > and R c 2 < 0, respectively (a), M = 1,3, 5, 7, 15 and N = 16 for i?§ > (b), M = 15 with 
N = 16,24 for R% > (c). 



along the x- or y-axis, which are both theoretically equivalent, this small difference is 
due to the six-fold rotational symmetry of the mesh, which is invariant upon rotation 
by 60° but not by 90°. For N — 24, no difference is noticeable between the R x > 1 and 
R x < 1 cases. 

Subsequently, we search for the disk radius in the following area-conserving form 



M+l 



R 2 (</>) = 1 + [ R °m cos(to0) + R s m sin(m0)], 



(4.2) 



where R c m and R s m are unknown amplitudes of cosine and sine terms in the Fourier 
series expansion of R 2 {(j)). Due to the area conservation and the mass centre fixed at 
the origin, (4.2) does not contain m — and m = 1 terms. Moreover, owing to the 
rotational invariance, we can set i?| = 0, which fixes the orientation of the drop up to a 
rotation by 90° provided that R% 7^ 0- This leaves 2M — 1 unknown coefficients in (4.2) 
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Figure 8. Equilibrium shapes for symmetries m = 2 (a) and m = 3 (b) at various magnetic 

Bond numbers. 

for the minimisation of the associated energy (2.20). The number of azimuthal modes M 
is chosen to ensure the convergence of equilibrium shapes. In this case, the mesh of unit 
circle is deformed radially to fit the disk 

(r i ,^ i ) = (r?ii(0?) ) 0?) ) (4.3) 

where (rf ,0?) and (r i; (pi) are the polar coordinates of the mesh points for circular and 
deformed disks, respectively. 

For comparison with the case of ellipse considered above, we start with M = 1, which 
leaves only one coefficient, i?2, in (4.2) to be determined. As seen in figure 6, which shows 
i?2 versus Bm for three numerical resolutions and two perpendicular orientations of the 
drop determined by the sign of the radial deformation of the mesh (4.3) results in 
a reduced numerical accuracy of the critical value of Bm, which for N — 24 elements in 
the radial direction is about 6% lower than its theoretical value (3.14). There is also a 
small difference in the shape depending on whether the drop is squeezed along the y- or 
x-axis (see figure 7a). Figures 7(a) and (b) show that the shape changes very little as the 
number of azimuthal modes and that of the elements in radial direction reach M = 15 
and N = 24, respectively. In the following, we will be using these values unless stated 
otherwise. 

The equilibrium shapes found as the magnetic field is gradually increased are shown in 
figure 8a. At Bm c rts 13.6, which due to the numerical approximation is slightly below the 
theoretically predicted stability threshold (3.14), the drop turns noticeably elliptic and 
rapidly elongates with a further increase in Bm. For Bm > 15, the drop starts to tighten 
around the middle part. No equilibrium shapes of this type can be found for Bm > 25. 
This implies that the drop may split up into two as the narrowing of the middle reaches a 
certain critical value. The splitting of the drop is not captured by this numerical method, 
which breaks down as the neck between two parts of the drop becomes too thin. 

Alternatively, when the magnetic field is applied instantly with Bm « 20 to a drop 
with some initial m = 3 perturbation, equilibrium shapes with a three-fold rotational 
symmetry shown in 8b are obtained. As seen in figure 9, the associated energy of m = 
3 mode is higher than that of m = 2 mode, which is also possible at the same Bm. 
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Figure 9. The associated energy (2.20) versus Bm for the shapes of various rotational symme- 
tries. Eo = 2-7T + %Bm is the associated energy of circular disk. The dots show the analytical 
bifurcation points (3.13). 
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Figure 10. Equilibrium shapes for symmetries m = 4 (a) and m = 5 (b) at various magnetic 

Bond numbers. 



Nevertheless, the shapes with three-fold symmetry are stable because they are separated 
from the two-fold symmetry shapes by a finite energy barrier. 

This, however, is not the case for the m — 4 and m = 5 symmetry shapes shown 
in figure 10(a) and (b), which can be obtained only when the corresponding symmetry 
is explicitly imposed in series (4.2) by ignoring all other modes. As seen in figure 11, 
the associated energy of four- and fivefold symmetries, in contrast to that of two- and 
threefold symmetries, decreases upon m = 2 and m = 3 radius perturbations. This 
implies that four- and fivefold symmetry shapes are indeed unstable with respect to 
these perturbations. 
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Figure 11. The associated energy normalised with its minimum for near equilibrium shapes 
of rotational symmetries with m = 2, 3, 4, 5 versus the perturbation of the cosine amplitude of 
m = 2 (a) and m = 3 (b) modes defined by R% and R% coefficients in (4.2). 



5. Summary and conclusions 

In this study, we have numerically modelled strongly deformed equilibrium shapes of 
a flat liquid metal drop subject to a transverse high-frequency AC magnetic field. The 
drop was treated as a thin liquid layer confined in a horizontal gap between two parallel 
insulating plates. AC frequency was assumed high so that the magnetic field was effec- 
tively expelled from the drop by the skin effect. Equilibrium shapes of the drop were 
found by using a variational principle for the associated energy involving the surface and 
magnetic contributions. Using Biot-Savart law, the associated electromagnetic problem 
was formulated in terms of a surface integral equation for the scalar magnetic poten- 
tial. This equation was solved numerically on an unstructured triangular mesh covering 
the surface of the drop. Numerical method was validated against analytical solution for 
the stability of circular disk with respect to small-amplitude azimuthally harmonic edge 
perturbations. According to the analytical solution, the edge deformations with the az- 
imuthal wavenumbers m = 2, 3, 4, . . . start to grow on circular disk as the magnetic Bond 
number exceeds the critical threshold values Bm^ = 3n(m+l)/2. The most unstable is 
m = 2 mode, which corresponds to an elliptical deformation at the critical Bond number 
5m< 2) = 9tt/2 « 14.1. 

This result agrees surprisingly well with the experimental findings of Conrath et al. 
(2006); Conrath (2007) for a drop of Galinstan (GalnSn eutectic alloy) with the diameter 
of 2R = 65 mm confined in a horizontal gap between two parallel glass plates separated 
by h — 3 mm. The drop was submitted to the AC magnetic field generated by a 10- 
winding (n = 10) coil with the inner and outer radii of R\ = 48 mm and i?2 = 81 mm, 
respectively, which were roughly in the plane of the drop. For the AC frequency of / w 
43 kHz, which was the highest one applied in the experiment, an originally circular disk 
became elliptical as the effective current in the coil exceeded / w 75 A. This corresponds 
to the r.m.s. magnetic field in the centre of the coil 

Bo Honl In R2 — In R\ 



V2 2 R2-R1 



7.5 mT, 



BR 

which yields the critical Bond number Bm = ~ 14, where 7 = 0.718 N/m is 

the surface tension of Galinstan and the effective arclength of the edge lo w nh/2 is 
approximated by a half circle. The critical currents are higher at lower frequencies and 
appear to saturate as the frequency is increased, which is consistent with the saturation 
of the electromagnetic force in the perfect conductor limit. Shapes with a rough three- fold 
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rotational symmetry are observed above the critical current I w 100 A, which corresponds 
to Bm w 25. This is by about a third greater than the theoretical value Bmf = 67r w 19 
for m = 3 mode. Note that also the shapes with a four-fold rotational symmetry are 
observed in the experiment though the numerical simulation showed them to be unstable. 
These discrepancies between the theory and experiment may be due to two effects. First, 
the size of the drop is comparable to that of the coil, which makes the applied magnetic 
field non-uniform over the drop radius. Second, to prevent the oxidation the drop is 
submerged in a 6% solution of HC1, which may affect the surface tension. Given all 
these experimental uncertainties and deviations from the idealised theoretical model, the 
agreement of the instability threshold for the m = 2 mode seems too good and perhaps 
even incidental. 

Note that the critical Bond number resulting from the energy variation approach is by 
a factor of 3 greater than that supplied by our previous linear stability analysis (Priede 
et al. 2006). There seem to be no obvious errors in either approach except for the factor 
of 2 missed in the final expression for the time-averaged force Fq above equation (24) of 
Priede et al. (2006). This factor taken into account results in Bm^ = 7r(m + l)/2 which 
increases the actual difference from 1.5 to 3 times. The only questionable point is the 
determination of electromagnetic force on the edge, where the magnetic field becomes 
singular, by the integration of Maxwell stress tensor over a small cylindrical surface 
enclosing the edge (Priede et al. 2006). It is important to notice that the local magnetic 
field at the edge used in the integration is entirely due to the currents induced in the 
sheet. Using Ampere's force law, it can be shown that such an approach accounts only 
for the interaction between the induced currents while it misses out any interaction of the 
induced and external currents. This is because the latter act via the external magnetic 
field, which is opposite to the induced one, but not taken into account by the local field 
distribution. As a result, the force on the edge is overestimated and, consequently, the 
magnetic field strength necessary for the instability underestimated. Obviously, the semi- 
infinite sheet model used by Priede et al. (2006) is not able in principle to account for the 
interaction with external magnetic field, which requires the consideration of finite size 
system. This is implied also by the energy variation approach, which does not work for a 
semi- infinite sheet model. On the one hand, the energy of the magnetic field, which falls off 
as ~ 1/\A* from the edge, diverges for semi-infinite sheet. On the other hand, this energy 
does not vary with the variation of the edge position because this variation is equivalent 
to the offset of the origin of coordinate system. This makes the force on the edge of semi- 
infinite sheet undetermined. Such ambiguities do not arise when the energy variation 
approach is applied to finite-size drops, as done in this study. Moreover, difficulties due 
to the edge singularity disappear altogether when smooth drops are considered, which, 
however, significantly increases the numerical complexity of the problem. 

I would like to thank Yves Fautrelle for stimulating discussions. 



Appendix A. The magnetic potential for harmonically deformed disk 

In the oblate spheroidal coordinates, (2.16) for the azimuthal mode m of the potential 
defined by (3.3) takes the form 
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The potential perturbation, which is supposed to vanish with the distance from the disk 
\&2™ — > 0, is related with the radius perturbation by (3.2), which now reads as 



*5" 



2 

TTT] 



(A 2) 



Although (A 1) admits the variable separation, such a solution is complicated by the edge 
singularity (A 2). Nevertheless, a compact analytical solution can be found similarly to the 
construction of spherical solid harmonics from the fundamental solution of the Laplace 
equation (Batchelor 1973) as follows. Firstly, note that if "J is a solution of the Laplace 
equation and e is a constant vector, then (e- V)\& is a solution, too. Secondly, if ^ satisfies 
a homogeneous boundary condition and e is directed along the boundary, then (e • V)^ 
satisfies that boundary condition, too. Thirdly, the operator (e • V) changes the radial 
dependence of 1 J> from <~ (r — 1)" to ~ (r — 1) Q_1 , while the azimuthal dependence is 
changed from the mode mtom+1. Algebra becomes particularly simple when e is taken 
in the complex form as e = e x + \e y = c 1 *(e r + ie^). Then each application of (e • V) 
is accompanied by the multiplication with e 1 ^. Thus, the solution for m = 1 is obtained 
straightforwardly from the axisymmetric base state (3.1) as 



*i(^) = -c- i0 (e-V)* o 



/I 2\ 1/2 

7T / 1 — 1} \ Tj 



2\i+e J n 2 +e' 



Higher azimuthal modes can be obtained similarly as = e vm ^ (e • V)" 1 ^f™, where 
&™ is an axisymmetric solution satisfying (Al). The edge condition (A 2) 

*™ 1 

n 2m 1] 

yields 'J'™ <~ rj 2 " 1 ^ 1 for 77 — ^ 0. Moreover, the perturbation vanishes far away from the 
disk when *I>™ = c^n 2 " 1 ^ 1 along the whole disk, where c™ is a constant. Then the 
corresponding axisymmetric solution of (Al) can be written as 



fc=l 



where P n (x) and Q n {x) are the Legendre polynomials and functions of the second kind, 
respectively (Abramowitz & Stegun 1972); the expansion coefficients are found as c™ = 

CT) 7 ™ where 



Jo 



r> 0F2 1 - 2m (2m-l)! 
P2k-i{V) d?7 = 



(m-k)W(m + k + 1/2) 
Then the solution for the perturbation amplitude can be written as 

*? = D+_ 1 D+_ 2 ---DtD a>V, (A3) 



using the operator 



n± = r ( F— - 1L\ + ™ 



which is defined by D± = e -^ m±1 ^ ( e± . V)e im<t >, where e+ = e, and e_ = e* is the 
complex conjugate of e. The calculation of (A3) is algebraically complicated but can 
be done by the computer algebra system Mathematica (Wolfram 1996), which requires 
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considerable computer resources and practically can be carried out only for m < 5. But 
this suffices to deduce the general solution (3.4). 

The axisymmetric solution (3.6) with ~ r/~ 3 edge singularity can be obtained in a simi- 
lar way directly from the axisymmetric base solution (3.1) by applying (e_ • V) (e+ • V) = 
D^Dq. This operator is equivalent to —dl because it represents the transversal part of 
the Laplace operator while (3.1) satisfies the Laplace equation. 
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